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Abstract. The reconstruction of the equilibrium of a plasma in a Tokamak is a 
free boundary problem described by the Grad-Shafranov equation in axisymmetric 
configuration. The right-hand side of this equation is a nonlinear source, which 
represents the toroidal component of the plasma current density. This paper deals 
with the identification of this nonlinearity source from experimental measurements in 
real time. The proposed method is based on a fixed point algorithm, a finite clement 
resolution, a reduced basis method and a least-square optimization formulation. This 
is implemented in a software called Equinox with which several numerical experiments 
are conducted to explore the identification problem. It is shown that the identification 
of the current density averaged over the magnetic surfaces is very robust. 
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1. Introduction 

In fusion experiments a magnetic field is used to confine a plasma in the toroidal vacuum 
vessel of a Tokamak p]. The magnetic field is produced by external coils surrounding 
the vacuum vessel and also by a current circulating in the plasma itself. The resulting 
magnetic field is helicoidal. 

Let us denote by j the current density in the plasma, by B the magnetic field and 
by p the kinetic pressure. The momentum equation for the plasma is 

du _ . 

9 ~dt + = J X 

where u represents the mean velocity of particles and p the mass density. At the slow 
resistive diffusion time scale [2J the term p— can be neglected compared to Vp and the 
equilibrium equation for the plasma simplifies to 

j x B = Vp 

meaning that at each instant in time the plasma is at equilibrium and the Lorentz 
force j x B balances the force Vp due to kinetic pressure. Taking into account the 
magnetostatic Maxwell equations which are satisfied in the whole space (including the 
plasma) the equilibrium of the plasma in presence of a magnetic field is described by 

p j = V x B, (1) 
V-B = 0, (2) 
j x B = Vp, (3) 

where po is the magnetic permeability of the vaccum. Ampere's theorem is expressed 
by Eq. (0Q) and Eq. (T5]) represents the conservation of magnetic induction. From the 
equilibrium equation (j3J) it is clear that 

B • Vp = and j • Vp = 0. 

Therefore field lines and current lines lie on isobaric surfaces. These isosurfaces form a 
family of nested tori called magnetic surfaces which enable to define the magnetic axis 
and the plasma boundary. On the one hand the innermost magnetic surface degenerates 
into a closed curve and is called magnetic axis and on the other hand the plasma 
boundary corresponds to the surface in contact with a limiter or to a magnetic separatrix 
(hyperbolic line with an X-point). 

The Grad-Shafranov equation [31 EJ E] is a rewriting of Eqs. (JIH) under the 
axisymmetric assumption. Consider the cylindrical coordinate system (e r , e^, e.,). The 
magnetic field B is supposed to be independent of the toroidal angle <p. Let us decompose 
it in a poloidal field B p = B r e r + B z e z and a toroidal field B^ = B^e^ (see Fig. [p. 

Let us also introduce the poloidal flux 



ib(r, z) = — / Bds = / B z rdr 



Figure 1. Toroidal geometry. 



where D is the disc having as circumference the circle centered on the Oz axis and 
passing through a point (r, z) in a poloidal section. From Eq. (T5]) one deduces that 

Bp = -[Vip x e^]. Therefore B.V^ = meaning that if) is a constant on each magnetic 

surface and that p = p(ip). 

The same poloidal-toroidal decomposition can be applied to j. From Eq. <^ it is 

clear that V • j = 0. As for B p it is shown that there exists a function /, called the 

1 / 

diamagnetic function, such that j„ = -[V( — ) x eJ. Since j.Vp = then V/ x Vp = 

r n G 

and / is constant on the magnetic surfaces, / = f{ip). 

f 

From Eq. ([I]) one also deduces that B^, = — and 
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and since 



Vp = p'(^)V^ and V/ = f(ip)Vip 
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the Grad-Shafranov equation valid in the plasma reads 

-AV = rp / (V') + — (ff')W (4) 

Thus under the axisymmetric assumption, the three dimensional equilibrium Eqs. 
(P - [3]) reduce to a two dimensional non linear problem. Note that the right-hand side 
of Eq. ([!]) represents the toroidal component of the current density in the plasma 
which is determined by the unknown functions p' and //'. In the vaccum there is no 
current and the poloidal flux statisfies 

-A*<0 = 

In this paper, we are interested in the numerical reconstruction of the equilibrium 
i.e of the poloidal flux ijj and in the identification of the unknown plasma current density 
[6l El Ej - In a control perspective this reconstruction has to be achieved in real time from 
experimental measurements. The main difficulty consists in identifying the functions 
p' and //' in the non linear right-hand side source term in Eq. (T4j). An iterative 
strategy involving a finite element method for the resolution of the direct problem and 
a least square optimisation procedure for the identification of the non linearity using a 
decomposition basis is proposed. 

Section 2 is devoted to the statement of the mathematical problem and to the 
description of the experimental measurements avalaible. The proposed algorithm is 
described in Section 3. This methodology has been implemented in a software called 
Equinox and numerical results using synthetic and real measurements are presented in 
Section 4. 



2. Setting of the direct and inverse problems 

2.1. Experimental measurements 

Although the unknown functions p'{ip) and (ff')(ip) cannot be directly measured in a 
Tokamak several measurements are available: 

• Magnetic measurements: they represent the basic information on which any 
equilibrium reconstruction relies. Flux loops provide measurements of ip and 
magnetic probes provide measurements of the poloidal field Bp at several points 
around the vacuum vessel. Let Q be the domain representing the vacuum vessel 
and 8Q its boundary. In what follows we assume that we are able to obtain the 

Dirichlet boundary conditions ip = 9d and the Neumann boundary conditions 

ldijj 

-—— = gN at any points of the contour dQ thanks to a preprocessing of the 
r on 

magnetic measurements. This preprocessing can either be a simple interpolation 
between real measurements or be the result of some boundary reconstruction 
algorithm which computes ip outside the plasma satisfying A*^ = under the 
constraint of the measurements [9j [10l [11] . 
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A second set of measurements which can be used as a complement to magnetic 
measurements are internal measurements: 

• Interferometric measurements: they give the values of the integrals along a family 
of chords Cj of the electronic density n e (ip) which is approximately constant on each 

flux line / n e (ip) dl = 7$. 

• Polarimetric measurements: they give the value of the integrals 

J Ci r on 

dip 

— — is the normal derivative of ijj along the chord q. 
on 

Even when using magnetic measurements only for the equilibrium reconstruction 
the numerical algorithm presented in this paper also uses: 

• Current measurement: it gives the value of the total plasma current I p defined by 

Ip = / j<pdx. 



Ampere's theorem shows that this quantity can be deduced from magnetic 
measurements. 

Toroidal field measurement: it gives the value B Q of the toroidal component of the 
field in the vacuum at the point (ro, 0) where ro is the major radius of the Tokamak. 
This is used for the integration of //' into / and for the computation of the safety 



factor q (see | Appendix A[). 



2.2. Direct problem 

The equilibrium of a plasma in a Tokamak is a free boundary problem. The plasma 
boundary is determined either as being the last flux line in a limiter L or as being a 
magnetic separatrix with an X-point (hyperbolic point). The region Q p C Q containing 
the plasma is defined by 

il p = {x e f2, ^( x ) > 4>b} 

where ipi, = maxLip in the limiter configuration or ip^ = ij}(X) when an X-point exists. 

In the vacuum region, the right-hand side of Eq. |4] vanishes and the equilibrium 
equation reads 

A*^ = in Q \ Q p 
_ — 

Let us introduce the normalized flux if) = — e [0, 1] in Q„ with 

ip a = maxfi 1/;, AN)) = -^-p'(ip) and B(ip) = (ff')Nj). This is introduced so 

A A/i r 
that the non dimensional and unknown functions A and B are defined and identified on 
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the fixed interval [0,1]. Imposing Dirichlet boundary conditions the final equilibrium 
equation is expressed as the boundary value problem: 

r — Tn — 

-A*V> = A[-A(V) + -B&)] X n P mfi 

r r ( 5) 

ip = Qd on dfl 

The free boundary aspect of the problem reduces to the particular non linearity 
appearing through xn p the characteristic function of Q p . The parameter A is a scaling 
factor used to ensure that the given total current value I p is satisfied 

I p = X [ [—A$) + ^B$)]dx. (6) 
Jn p r r 

2.3. Inverse problem 

The inverse problem consists in the identification of functions A and B from the 
measurements available. It is formulated as a least-square minimization problem 

Find A*, B*, n* e such that : 

(7) 

J {A*, B\ n*) = inf J(A, B, n e ). 

If magnetic measurements only are used the formulation only needs the A and B 
variables and the J\ and J<i terms in Eq. (jSj) below are not needed. When polarimetric 
and interferometric measurements are used, the electronic density n e (ip) also has to be 
identified even if it does not appear in Eq. (jSJ). The cost function J is defined by 

J (A, B, n e ) = J + J 1 + J 2 + J £ (8) 
Jo describes the misfit between computed and measured tangential component of B p 

k=l 

where N is the number of points of the boundary dfl where the magnetic 
measurements are given. 

N, 



T 1 \~~V Volar\2i f n e(V ; ) dtp , 2 



k = l JC k 

and 

N, 



J2 = lY,H nter ) 2 (f n e $)dl-lk) 



k=i " Ck 

N c is the number of chords over which interferometry and polarimetry 
measurements are given. The weights w give the relative importance of the different 
measurements used. As a consequence of the ill-posedness of the identification of A, B 
and n e , a Tikhonov regularization term J e is introduced [12] where 

Js = y j\A"{x)fdx + £ -f J\b"{x)?oIx + £ -f j\K{x)fdx 
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and Eai £b and e ne are the regularization parameters. 

The values of the different weights and parameters introduced in the cost function 
are discussed in Section HJ 

It should be noticed here that magnetic measurements provide Dirichlet and 
Neumann boundary conditions. The choice was made to use the Dirichlet boundary 
conditions in the resolution of direct problem and to include the Neumann boundary 
conditions in the cost function formulated to solve the inverse problem. This is arbitrary 
and another solution could have been chosen. 



3. Algorithm and numerical resolution 

3.1. Overview of the algorithm 

The aim of the method is to reconstruct the equilibrium and the toroidal current density 
in real time. At each time step determined by the availability of new measurements 
during a discharge, the algorithm consists in constructing a sequence (ip n , flp, A n , B n , A n ) 
converging to the solution vector (ip, Q p , A, B, A). The unknown function n e may be 
added too if interferometry and polarimetry measurements are used. The sequence is 
obtained through the following iterative loop: 

• Starting guess: Q®, A , B° and A known from the previous time step solution. 

• Step 1 - Optimisation step: compute A n+1 satisfying ([6]) 

A™ +1 = IJ [ {^-A n (4> n ) + —B n (tp n )]dx 

then compute A n+1 (ip n ) and B n+1 (ijj n ) using the least square procedure detailed in 
Section 13X21 

• Step 2 - Direct problem step: compute of i/j n+l solution to 
_ A *^n+i = X n+1 [^A n+ \$ n ) + ^B n+1 (i) n ))xn$ infi g 

tp n+1 = 9d on d£l. 
and the new plasma domain 

• n :— n+ 1. If the process has not converged return to Step 1 else (ip, fl p , A, B, A) = 
(ip n , flp, A n , B n , A n ). The process is supposed to have converged when the relative 

residu r . — r . is small enough. 

Il^ll 

At each iteration of the algorithm, an inverse problem corresponding to the 
optimization step and an approximated direct Grad-Shafranov problem have to be solved 
successively. In Eq. (Q, tp n is known and since the right-hand side does not depend on 
ip n+1 the boundary value problem is linear. 

In the next section the numerical methods used to solve the two problems 
corresponding to step 1 and step 2 are detailed. 
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3.2. Numerical resolution 

3.2.1. The finite element method for the direct problem The resolution of the direct 
problem is based on a classical P l finite element method [13]. Let us consider the family 
of triangulation of Q, and Vh the finite dimensional subspace of H 1 ^) defined by 

V h = {v h G H\n),v hVr E P\T), VT G r h }. 

and introduce V h ° = 14 fl Hq(Q). The discrete variational formulation of the boundary 
value problem [9] reads 

Find iph G Vh with ip h = g D on <9f2 such that 

r 1 /■ r _ k> (10) 

Vv fc G / ■ Vv h dx = / A[— A($*) + ^(V^W^ 

/^o^ Jn p Mo r 

where ip* represents the known value of ip at the previous iteration. Numerically the 
Dirichlet boundary conditions are imposed using the method consisting in computing 
the stiffness matrix K of the Neumann problem and modifying it. Consider (vt) a 

f 1 

basis of Vh then K t j = / VviVvjdx The modifications consist in replacing the rows 

Jn t l o r 

corresponding to each boundary node setting 1 on the diagonal terms and elsewhere. 
At each iteration only the right-hand side of the linear system in which the Dirichlet 
boundary conditions appear has to be modified. The linear system corresponding to 
Eq. [TU] can be written in the form 

K.* = y + g (11) 

where K is the n x n modified stiffness matrix, \I/ is the unknown vector of size 
n the number of nodes of the finite elements mesh, y is the vector associated with the 
modified right-hand side of Eq. fflOl) and g is the vector corresponding to the Dirichlet 
boundary conditions. 

The matrix K is sparse and so is its LU decomposition. The inverse matrix K~ x 
however is not sparse. The linear system flllD is inverted using the LU decomposition 
since it is computationally cheaper than using the full inverse matrix K" 1 which is 
nevertheless needed for the optimization step of the algorithm in Eq. (JTSJ) below. 

The vector y depends on functions A and B which are determined in the 
optimization step. Functions A, B and n e are decomposed on a finite dimensional 
basis (3>i)i=i,..., m of functions defined on [0, 1] 

mm m 

A(x) = ^^aj$j(x), B(x) = ^^6,$i(x) and n e (x) = Cj$j(x). 



The vector y reads 

y = Y{%jj*)u 



(12) 
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where u = (a±, ...,a m ,bi, ...,b m ) G M? m is the vector of the components of functions A 
and B in the basis The matrix Y of size n x 2m is defined as follows. Each row % 
of Y is decomposed as 



3.2.2. Detailed numerical algorithm One equilibrium computation corresponds to one 
instant in time during a pulse. The quasi-static approximation consists in considering 
that at each instant the Grad-Shafranov equation is satisfied. During a pulse successive 
equilibrium configurations are computed with a time resolution At corresponding to the 
acquisition time of measurements: 

• Initialization before the discharge: the modified stiffness matrix K, its LU 
decomposition as well its inverse K~ x are computed once for all and stored. 

• Consider that the equilibrium at time t — At is known and that a new set of 
measurements is acquired at time t. 

• Computation of the new equilibrium at time t through the iterative loop briefly 
described in the previous Section and detailed below: 

The equilibrium from the previous time step is used as a first guess in the iterative 
loop. 

Step 1 - Optimization step During the optimisation step, n e is first estimated from 
interferometric measurements and A and B are computed in a second time. 

• Compute the electronic density n e based on the equilibrium of the previous 
iteration ip* using a least square formulation for the minimun of Ji with Tikhonov 
regularization and solving the associated normal equation: The flux ijj* is given. 




TO 




The interferometric measurements for i — 1 ... n c are 




The cost functional reads 



J{v) = ^YM^fCZ.B^-^ + ^i^v 



= lllD^Bv-^ + ^Av 
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where D 1 ! 2 = diag(wl nter ) and the regularization matrix A is defined by 

Jo 

and $• is the second derivatives of the basis function 
It is minimized solving the associated normal equation 

(a\D l ' 2 B) T {D l ' 2 B) + eA)v = a(D 1 ' 2 B) T D 1 '^ (13) 

Since n e ~ 10 19 m -3 an adimensionalizing parameter a = 10 19 m -3 , such that 
v = av, is introduced in order to precondition the linear system which is inverted 
using LU decomposition, as well as a reasonable prescribed value for the non 
dimensional regularization parameter e = a 2 e. 

Compute A n+1 satisfying Eq. ([6]). In the right-hand side y, A appears in the product 

Am. In order to avoid any divergence issue due to the non uniqueness of A (for all 
u 

a, Xu = (Aa)(— )) the degrees of freedom (dofs) u are scaled by m = max(laJ), u 
1 

is replaced by — u and A by mX. 

m 

Compute A and B. In order to approximate A and B, suppose n e is known and 
consider the discrete approximated inverse problem 

Find u minimizing : 

j( U ) = 1 -\\c(r)*-d\\i + £ -u T Au (14) 

where C(ip*) is the observation operator and d the vector of experimental 
measurements. The first term in J is the discrete version of Jo + J\. The second 
one corresponds to the first two terms of the Tikhonov regularization J £ with 
sa = £b = £ which will always be assumed in order for functions A and B to 
play a symmetric role. 

Let us denote by I the number of measurements available (I — N + N c if magnetic 
and polarimetric measurements are used) and by D the diagonal matrix made of 
the weights and wl° lar , the norm is defined by Vx £ Mf ||x|||, = (Dx,x) = 

C{ip*) is a sparse matrix of size / x n and can be viewed as a vector composed 
of two blocks Co of size N x n and indepedent of ip* and Ci(ip*) of size N c x n 
corresponding respectively to J and J\. That is to say that multiplication of the 
kth row of Co by ip gives the kth Neumann boundary condition approximation 

1 dtb 

Co *¥ « (-IT (M* . 
r on 

The block C\(ip*) depends on ip* through the n e (ip*) function. The multiplication of 
the kth row of Ci(ip*) by ^ gives the kth polarimetric measurements approximation 
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The matrix A is of size 2m x 2m and is block diagonal composed of two blocks Ai 
and A 2 of size m x m, with 

(Ai)y = (A 2 )y = / <I>- / (./ )<I>- / (./)r/./ 
./0 



Using Eqs . (fill - [121) problem ffl4|) becomes 

= ii^OTir^c^ti + {CW)K- X g - d)f D + £ -u T Au 

= l\\Eu-f\\l+ £ -u T Au 

where E = C^K^Y^*) and / = -C^K^h + d. Setting E = D l ' 2 E, 
problem tHM reduces to solve the normal equation 

(E T E + eA)u = E T f (15) 

whose solution is denoted by u*. 

Direct problem step Update the dofs u and update the flux ip by solving the linear 
system 

Kip = Y(4)*)u* + g 

using the precomputed LU decomposition of matrix K. Update fl p possibly computing 
the position of the X-point if the plasma is not in a limiter configuration. 

4. Numerical results 

4-1. Twin experiment with noise free magnetic measurements 

In this section we assume that the poloidal flux corresponding to an equilibrium 
configuration ip is given on the boundary V . These Dirichlet boundary conditions 
can either be real measurements or can be the output from some equilibrium simulation 

code. In a first step we also assume to know functions p 1 and //' (or A and B). It is 

dip 

then possible to run a direct simulation to compute tp on Q (see Fig. [2]) and thus ^r— 

ran 

on T which can then be used as measurements in an inverse problem resolution. 

In this first experiment the magnetic measurements are free of noise. The first guess 
unknown functions are A(x) = B(x) = 1 — x and A = 1. The poloidal flux ip is initially a 
constant on Q. The weights in the misfit part of the cost function J related to magnetic 

measurements are defined by Wk = . ■ Since the error on magnetic measurements 

are of about one percent we define a = 0.01B m where B m is a mean magnetic field value 



which thanks to Ampere's theorem can be defined as B n 
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Figure 2. An equilibrium configuration for the tokamak JET from which twin 
experiments are performed. The domain f2 is shown. Isofiux are plotted from ip = 
(magnetic axis) to ip = 1 (plasma boundary defined by the existence of an X-point at 
point r — 2.5 and z — —1.4 m) by step of Aip = 0.1 Interferometry and polarimetry 
chords appear in green. 



The functions A and B are decomposed in a function basis denned on the interval 
[0,1]. Several basis have been tested (piecewise affine functions, polynomials, Bplines 
and wavelets) in order to verify that the result of the identification does not depend on 
the decomposition basis. This is the long as the dimension of the basis is large 

enough. In the remaining part of this paper each function is decomposed in the same 
basis of 8 Bsplines [14] . The boundary condition A(l) = B(l) = is imposed. 

The computations are carried out for several values of the regularization parameters 
e ranging from 10~ 10 to 1. We are interested in the ability of the method to recover 
functions A and B and thus the current density profile averaged over the magnetic 



surfaces (see Appendix Appendix A ): 



ro < >= + Ar 2 < 1 > B$) 



and the safety factor q (see | Appendix B[ ). 



As can be seen from Fig. [3] the optimal choice for e is of about 10 5 for which 



13 



functions A and B are well recovered. For smaller values some oscillations appear 
because the regularization is not strong enough and on the contrary greater values lead 
to less precision in the recovery of the unknown functions since regularization is too 
strong. In the second column the relative errors on the identified functions are plotted. 




Figure 3. Twin experiment with noise free measurements and different regularization 
parameters e ranging from 1CP 10 to 1. Left column: identified functions XA{^i), 

Atq < — > B( , 0) for each different e value, and the known reference functions. Right 
column: relative errors. 

Figure H] shows an important point. Almost whatever the chosen value of e is, i.e. 

whatever the quality of the identification of A and B is, the identified average current 
j(r, ip) 

density r < > as well as the safety factor q are always well recovered and the 

r 

relative errors are one order of magnitude smaller than for functions A and B. The same 
kind of observation was made in [8] where the identified functions A and B seemed to 
be rather sensitive to perturbations whereas the mean current density was very stable. 

In Table (TJ the evolution of the relative residu on ip, A, B and A versus the number 
of iterations is given. It demonstrates numerically the convergence of the algorithm 
in this case where a value of 10~ 6 is used as stop condition. The algorithm needs 10 
iterations to converge. It is interesting to notice that even though the first guess is 
not particularly well chosen the relative residu on ip at the second iteration has already 
fallen to 4%. In real applications when simulating a whole pulse the first guess for the 
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Figure 4. Twin experiment with noise free measurements and different regularization 

parameters e ranging from 10 -10 to 1. Left column: resulting identified mean current 
j (f^ ipj 

density ro < ■ — >, safety factor q for each e value and the corresponding known 

r 

reference values. Right column: relative errors. 



computation of the equilibrium at t is the equilibrium computed at t — 5t and 2 iterations 
are enough to ensure a good convergence of the algorithm. 

4-2. Twin experiment with noisy magnetic measurements 

Figures [5] and [6] shows the results of the same type of numerical experiment but with 

_ , , . dxj) 

noisy measurements. Each magnetic input, m representing either ip or — — at a point 

ran 

of the domain boundary T is perturbated with a one percent noise normally distributed, 
m v = m + rj with rj ~ N(m, 0.01m). For each chosen value of the regularization 
parameter the algorithm is run 200 times with measurements randomly perturbated as 

above. Then for each function XA, Xr^ < — > B, ro < — > and q, a mean 
function, a median function and a standard deviation function is computed. 

In comparison with the noise free case the regularization parameter needs to be 
significantly increased to values of at least e = 10~ 2 and for a safer convergence of the 
algorithm to e — 10 _1 . For smaller values the algorithm either does not converge or 
gives very oscillating identified functions. 
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Table 1. Numerical convergence of the algorithm. 



nTT 



nTT 



A n+i - A T 



Iteration n 



A 



A'' 



B 



A< 



A' 



1 


9 (SzLSDQ 


U.U t OJJ 


O.GOUV 


n 1 nm 97 

U . 1UU 1 


2 


0.0408642 


1.19473 


1.42619 


9.24968 


3 


0.0733385 


1.83005 


1.47338 


0.563235 


4 


0.0404254 


0.884617 


1.0359 


0.108107 


5 


0.00539736 


4.79091 


4.37571 


0.826455 


6 


0.000349811 


0.127626 


0.180449 


0.0889022 


7 


1.58606e-05 


0.0262942 


0.0246657 


0.0263 


8 


5.67036e-06 


0.00294791 


0.0024952 


0.00315952 


9 


1.4533e-06 


0.000339986 


0.000273055 


0.000362224 


10 


6.19066e-07 


6.41923e-05 


6.51076e-05 


6.29838e-05 



The mean error on the reconstructed functions is always smaller in the interval 
ip G [0.5,1] than in the interval [0,0.5]. This is due to the fact that magnetic 
measurements are external to the plasma and do not provide enough information to 
properly reconstruct the functions in the innermost part of the plasma. 

As e increases the variability or the standard deviation on the identified functions 
decreases. With small e the algorithm can find very different functions depending on 
the perturbations of the measurements. With e = 10~ 2 the variability in the identified 
functions A and B is strong however the mean identified functions are close to the exact 
reference ones. On the other hand with e = 1 the variability of the identified functions 
is strongly reduced but they are quite different from the exact reference functions in the 
interval [0, 0.5]. 

It is worth noticing that in all cases the resulting safety factor q and mean current 

j (r, tp) 

density r$ < — - — > are well recovered. The remark of the preceding section on 

r 

the identifiability of the mean current density still holds: it is quite well recovered 
even if functions A and B taken separately are not well identified. The mean error on 
the current density profile is almost always smaller than the mean errors on functions 
A and B. Moreover this error does not change very much between the different 
cases and particularly between the e = 10 _1 and the e = 1 cases. This implies 
that for a large interval of e the value of the part of the cost function related to 
magnetic measurements Jo is almost constant. Therefore it is difficult to find an optimal 
value for the regularization parameter. For example the L-curve method [15] for the 
determination of the regularization parameter can hardly be used and gives some results 
which are not very reliable since the L-curves are not well behaved and the location of 
the corner is not clear. The "L" is an almost vertical line. This due to the fact that 
in a large interval of e values an increase in e implies a important decrease in the 
regularization term ^(u*(e)) T Au*(e) but does not lead to a significative increase in the 
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misfit term J (u*(e)). 




median 



0.2 0.4 0.6 0.8 1 0.2 0.4 0.6 0.8 1 




0.2 0.4 0.6 0.8 1 0.2 0.4 0.6 0.8 1 



Figure 5. Statistical results of the identification experiments with noisy magnetic 
measurements. Row 1: e = 10~ 2 , row 2: e — 10 _1 , row 3 e = 1. Column 1: function 

XA(^j) and column 2: Xrfi < — > B(ip). For each function the reference value from 
which the unperturbated measurements were computed is given in black and the mean 
identified function in red. The mean ± standard deviation functions are shown in 
dashed red. 
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Figure 6. Statistical results of the identification experiments with noisy magnetic 
measurements. Row 1: e = 10 -2 , row 2: s — 10 _1 , row 3 e = 1. Column 1: 
j(v ib) 

ro < ■ — >, and column 2: safety factor q. For each function the reference value is 

r 

given in black and the mean identified function in red. The mean ± standard deviation 
functions are shown in dashed red. 
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4-3. Twin experiment with noisy magnetic, interferometric and polarimetric 
measurements 

In this last twin experiment interferometric and polarimetric measurements are also 
used. At first a density profile is prescribed, n e (x) = \f\ — x on [0,1], as well as the 
same A and B functions as in the previous twin experiments. Then similar to the 
preceding section the equilibrium is computed from given Dirichlet boundary condition. 
A set of artificial magnetic, interferometric and polarimetric measurements is generated. 
Finally several twin experiments with a 1% noise are performed and some statistics are 
computed. The weights related to interferometric and polarimetric measurements in the 
cost function are defined as 



w 



polar 
k 



1 



/Jf a polar 



-, with a polar = 10 1 radians 



• w inter = with a inter = 10 18 m" 3 

\fW c o- inier 

The determination of the regularization parameter for the density function n e is 
far less a problem than for functions A and B since for example the L-curve method 
works quite well in this case (see Fig. [10] in the next Section) and the n e function is well 
recovered as shown on Fig. [9j The regularization parameter for the density function is 
set to e ne = 10~ 2 . 

The statistical results of the twin experiments are shown on Figs. [7] and [8] for 3 
different values of e. The use of interferometric and polarimetric measurements adds 
supplementary constraints on the A and B functions. The variability in the recovered 
functions is less important than in the case where only magnetics are used particularly 
for -0 G [0, 0.5]. This is not surprising since the new measurements are internal and bring 
some information contained inside the plasma domain. Nevertheless it is not enough to 
perfectly reconstruct independently the A and B functions. This does not prevent an 
excellent recovery of the mean current density profile and of the safety factor q. This 
phenomenum already observed in the magnetics only case is emphasized here where the 
variability of the recovered profiles has decreased. 



4-4- A real pulse 

The algorithm detailed in this paper has been implemented in a C++ software called 
Equinox developed in collaboration with the Fusion Department at Cadarache for Tore 
Supra and JET. Equinox can be used on the one hand for precise studies in which the 
computing time is not a limiting factor and on the other hand in a real-time framework 
to reconstruct the successive plasma equilibrium configurations during a whole pulse. 
For the time being it is used on JET and ToreSupra pulses but can potentially be used 
on any Tokamak. 

During the real time analysis of a whole pulse an equilibrium is reconstructed 
from new measurements with a time step of At = 100 ms. For each equilibrium 
reconstruction the number of iterations of the algorithm is set to 2. This enables fast 
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Figure 7. Statistical results of the identification experiments with noisy 
measurements (magnetics, interferometry and polarimetry) . Row 1: e = 10~ 2 , row 2: 

e = 10 _1 , row 3 e = 1. Column f: function \A(ip), and column 2: Atq < — > Bfy). 
For each function the reference value from which the unperturbated measurements 
were computed is given in black and the mean identified function in red. The mean ± 
standard deviation functions arc shown in dashed red. 



enough computations while a very good precision is achieved since the initial guess for 
an equilibrium computation at time t is the equilibrium computed at time t — At. After 
1 iteration a typical value for the relative residu on ip is of 10 -2 and it is of 10 -3 after 2 
iterations. Table [2] gives the size of the finite elements mesh used at ToreSupra and at 
JET as well as typical computation times on a laptop computer. 

During the computations the expensive operations are the updates of matrices C 
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Figure 8. 



Statistical results of the identification experiments with noisy 



measurements (magnetics, interferometry and polarimetry) . Row f: s = 10 , row 



2: e = 10 , row 3 e = 1. Column 1: tq < 



>, and column 2: safety factor q. 



For each function the reference value is given in black and the mean identified function 
in red. The mean ± standard deviation functions arc shown in dashed red. 



and Y as well as the computation of products CK~ l and CK~ X Y . In order to reduce 
computation time the K~ x matrix is precomputed and only the ^-dependent part of 
C is delt with. The resolution of the direct problem, Eq. (fTTI) . is cheap since the LU 
decomposition of the K matrix is also precomputed. 

The choice of the regularization parameters is crucial since it determines the balance 
between the fit to the data and the regularity of the identified functions. It is also 
difficult as is shown in the preceeding section. Ideally they should be determined for 
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Figure 9. Statistical results for the identification of the denisty function n e with noisy 
intcrferometric measurmcnts. 





ToreSupra 


JET 


Finite element mesh 


Number of triangles 


1382 


2871 


Number of nodes 


722 


1470 


Computation time (1.80GHz) 


One equilibrium 


20 ms 


60 ms 



Table 2. Typical mesh size and computation time for ToreSupra and JET 

each equilibrium reconstruction. However this is not possible in a real-time application 
and the regularization parameters have to be set apriori to a constant value. From the 
twin experiments presented in the preceeding sections it is quite clear that a good value 
for e the regularization parameter is in the range [10~ 2 , 1]. By trial and error on different 
pulses at JET using magnetics, interferometry and polarimetry, it appeared that a value 
of e — 5.10 -2 gave good results. 

As for the identification of the A and B functions the choice of a good regularization 
parameter for n e is crucial. However in this case the L-curve method works quite well 
and it was used to determine the regularization parameters e ne a priori on a number of 
equilibriums for a few shots. The obtained values showed little variation and the choice 
of a mean value e = 0.01 proved to be efficient. Figure (TU] shows an example of an 
L-curve computed for the identification of n e . 

An example of the outputs from Equinox is presented in Figure [TTJ It is the 
equilibrium computed at time 44.5 for JET pulse number 77601 using either magnetic 
measurements only or both magnetic and polarimetric measurements. In both cases the 
density n e is computed from interferometric measurements. One can observe the position 
of the plasma in the vacuum vessel. Isoflux lines are displayed from the magnetic axis 
to the boundary. The interferometry and polarimetry chords are displayed. At JET the 
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L-curve 

3i 1 1 1 




Figure 10. Typical Lcurve for the determination of £ ne . It is a plot of the parametric 
curve x(e ne ) = log(^\\D^ 2 (Bv*(e n J - 7 )|| 2 ), y(e n J = log(±(v* (e n J) T Av* (e n J) 
where v*(e ns ) is the solution to Eq. (jT3|l . Hansen's algorithm [15] locates a corner at 
e„ e = 0.01. 

chords are ordered as follows: number 1 to 4 are the vertical chords from left to right, 
and number 5 to 8 are the horizontal ones from bottom to top. The default configuration 
which is used in this real-time case is to use chords 3 to 8 to identify the density n e and 
to use chords 3, 5 and 7 for polar imetry. The other chords are not allways reliable and 
might give poor quality measurements which can forbid the computation of a satisfying 
equilibrium. For each chord (used or not) the error between computed and measured 
interferometry is given in purple. These errors are about 1% for the active chords. The 
polarimetry absolute errors are given in yellow. Different graphs are plotted on the left 
of the display. On the first row the identified function A, and corresponding functions 
p' and p. On the second row the identified function B and corresponding function //'. 
The third row gives the toroidal current density in the equatorial plane and the fourth 
one shows the safety factor. Finally on the fifth row the identified n e function is plotted. 

One should notice that in the particular case shown here the introduction 
of polarimetry measurements enables the obtention of what is called a reverse q 
profile (i.e a non monotonic q profile) and the existence of which is confirmed 
by magnetohydrodynamic analyses (D. Mazon and JET contributors personal 
communications) . 

5. Conclusion 

We have presented an algorithm for the identification of the current density profile 
in the Grad-Shafranov equation and the equilibrium reconstruction from experimental 
measurements in real time. We have shown thanks to several twin experiments that even 
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though the unknown functions A and B (or p' and //') taken separately might not be 
allways exactly identified the resulting mean current density and safety factor seem to 
be allways well identified. We have also shown that the use of internal polarimetric 
measurements improves the quality of the identification but is still not enough to 
perfectly identify both A and B. Finally we have introduced the software Equinox 
in which this methodology is developed. 
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Figure 11. Graphical outputs from Equinox. Reconstructed equilibrium at time 44.5 
s for JET pulse number 77601. Top: magnetics only are used. Bottom: magnetics and 
polarimctry. See text for mor details. 
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Appendix A. Average over magnetic surfaces 

The method of averaging over the magnetic surfaces is detailed in [16] (p 242). The 
average < A > of an arbitrary quantity A on a magnetic surface S is defined as 



where V is the volume inside S. This notion of average has the following property: 



where is a closed contour ip = cte G (0, 1) and B p = 
Appendix B. Safety factor q 

The safety factor is so called because of the role it plays in determining stability (pQ p 
111). It can be seen as the ratio of the variation of the toroidal angle needed for one 
magnetic field line to perform one poloidal turn. 



Since q is the same for all magnetic field lines on a magnetic surface it is a function of 
ip (or ip). The expression of q used for computations is the following 




< A > 




1 





where C, 



is a closed contour ip = cte G (0, 1), — — and 
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